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Abstract 


The vibration suppression effectiveness of a flexible in-line marine cable vibration 
absorber is studied. The transfer matrix method is used to build various numerical models 
of vibration absorbers in marine cable systems. The models determine cable system 
natural frequencies, mode shapes and modal damping ratios. 

The introduction of absorber damping is shown to result in complex roots to the 
modal characteristic equations. A computer complex root solver is used to solve for the 
complex roots of the characteristic equations, resulting in complex system natural 
frequencies. The significance of complex natural frequencies is explained. Complex 
natural frequencies are used to calculate modal damping ratios. 

The models demonstrate that absorber effectiveness is heavily dependent on 
absorber location, absorber mass and absorber length. Parametric variation is used to 
achieve maximum effectiveness of the flexible in-line absorber. Even under optimum 
conditions, it is shown that the absorber provides insufficient damping to reduce vortex- 
induced vibrations in water. 

The same transfer matrix method ts used to evaluate the effectiveness of a mass- 
spring-dashpot type absorber in a marine cable system. This type of absorber is shown to 
produce adequate damping to reduce vortex-induced vibrations in water. The transfer 
matrix method used in this thesis is validated by analyzing the same system using an 
approach by Den Hartog [1]. The transfer matrix approach combined with complex root 
solving capability is shown to provide an effective analysis method for marine cable 
systems. 
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Chapter 1 


Introduction 


Long slender marine structures such as cables and pipelines are subject to a 
condition called vortex induced vibration. This occurs when the cylindrical structure 1s 


subjected to a cross flow current as shown in figure 1-1: 


Flow speed, U (m/s) o~ \ nf a 
__ 
— Ney, 
alan 
Vortices in wake 
Figure 1-1: Vortex Shedding 

At a critical value of flow speed, U, a periodic shedding of vortices occurs on the upper 
and lower surfaces of the cylinder. This behavior is known as a Von Karman street. As 
each vortex is shed into the wake, a localized pressure drop occurs. The result ts a 
harmonically varying vertical (lift) force on the cylinder. The frequency at which vortex 


shedding occurs 1s related to cylinder diameter and flow velocity by a dimensionless 


parameter called the Strouhal number: 


i 
St = U Gan) 
where: f = vortex shedding frequency (Hz) 


D = cylinder diameter (m) 


U = flow speed (m/sec) 





The Strouhal number is approximately 0.2 for a wide range of Reynolds numbers. 
Vortex-induced vibrations are especially severe when the vortex shedding frequency is 
near a resonance frequency of the structure. Cross-flow displacement of the cylinder at 
this “lock-in” frequency may be on the order of one cylinder diameter. Such vibrations 
can fatigue marine structures and cause eventual failure. 

This problem motivates the investigation into a dynamic absorber for long 
cylindrical marine structures. Before analyzing the effect of various damping devices, it is 
worth while to determine the magnitude of damping which must be introduced into a 
marine structure to cause significant reduction in vortex induced vibrations. In figure 1-2 
below, Griffin [13] shows how cross flow displacement varies with a dimensionless 
parameter called the reduced damping, Sg. Reduced damping 1s the ratio of damping 


force to fluid exciting force. 
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Figure 1-2: Vibration Displacement vs. Reduced Damping 








where: Se= Sa Se Gaiiiey OF D’)= reduced damping (1-2) 

S, = Strouhal number 

m, = linear density of cable (kg/m) 

Or= fluid density (kg/m’) 

D = cable diameter (m) 
The reason displacement amplitude approaches a plateau value at low reduced damping 
(values below Sg about 0.2) is due to the self-limiting nature of vortex induced vibrations. 
At values of displacement amplitude greater than approximately one cylinder diameter, 
the extreme motion of the cylinder begins to disturb the formation of the Von Karman 
street. Stuctural damping ts not the limiting factor at low values of reduced damping. The 
reduced damping 1s valuable in determining the amount of damping which must be added 
to the structure to have an effect on displacement amplitude. It would not be worthwhile 
to add damping such that reduced damping would not exceed the “knee in the curve” 
value of 0.2 illustrated in figure 1-2. This is also a very convenient approach since the 
damping ratio used to calculate reduced damping is the damping which one would 
measure in an in vacuo transient decay test. Thus, various absorber configurations can be 
evaluated without the need to consider hydrodynamic effects. 

This thesis will study various vibration absorber configurations for marine cables. 

Most of the analysis will focus on the type of absorber shown in figure 1-3 on the 


following page: 





damper 


A 


cable = ; 


en 


Rigid masses 


Figure 1-3: Cable Vibration Absorber 


This absorber was first proposed by Vandiver and Li [10]. Energy dissipation occurs by 
way of relative angular motion at the pivot between the two rigid masses. The damping 
device is integral to the pivot but it is shown schematically in figure 1-3. Butler [7] 
showed this type of damping device to be possible for marine cable applications. It was 
seen as a very versatile configuration for marine cables since such an absorber could be 
clamped on at a desired location while deploying a cable. Or, if the absorber were 
permanently installed, its size and orientation might permit the absorber to be reeled onto 
a cable drum when the cable 1s not deployed. 

Chapter 2 describes the transfer matrix approach for analyzing multi-component 
mechanical systems with finite boundaries. It is a very convenient method for studying 
the dynamic behavior of marine cable systems with attached damping devices. The 
analysis utilizes the absorber transfer matrix derived by Li [5]. 

In chapter 3, the transfer matrix method 1s used to determine natural frequencies, 
mode shapes and damping ratios of the cable / absorber system. A method 1s formulated 
for determining the optimum value of damping constant. 

In chapter 4 , system parameters such as absorber length and mass are varied to 
study their effect on damping performance. Also, use of multiple absorbers and the effect 


of absorber location relative to vibration nodes and antinodes was studied. Use of a more 





conventional absorber type 1s studied by merging the transfer matrix approach with a more 
classical approach (Den Hartog [1]) with some interesting results. 


Chapter 5 concludes the study and suggests areas for further research. 
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Chapter 2 


The Transfer Matrix Method 


The transfer matrix method allows for a convenient matrix representation of the 
marine cable system. In 1963, Pestel and Leckie [2] recognized that the introduction of 
numerical computing made matrix methods a powerful tool for studying mechanical 
system dynamics. In marine cable systems, features such as rigid bodies and discrete 
damping devices are well suited to matrix modeling techniques. Today’s computing speed 


allows complicated marine cable systems to be modeled with relative ease. 


2.1 The Marine Cable System 


2.1.1 Assumptions 
It 1s important here to state the assumptions employed in representing the cable 

system in this way. These assumptions are those of the classic linear string wave equation. 
They were used by Li [5] in developing the transfer matrices used in this analysis: 

1) the cable systems are taut and undergo only small deflections 

2) cable bending stiffness is neglected 

3) the masses in the cable system are rigid 

4) the weight of these masses is small compared to cable tension 

5) the dynamic variation in tension in the cable is small and neglected 


6) only planar motion is considered (2 dimensional) 


1] 





2.1.2 Transfer Matrix Representation of a Marine Cable System 
Transfer matrices are used to relate the physical state of one location in a 
mechanical system to that at another location in the system. The state at any location in 
the system 1s represented by the state vector. In a cable system, the state at any location 

can be represented by the state vector, z: 


y(x) 
= a 
Ty (x) Ge) 


[z] 


where y is the Fourier transform of cable transverse displacement and y” 1s the Fourier 
transform of cable slope. The product of cable tension T and y’ ts a value equal to the 
transverse force caused by tension at the location of the state vector. Thus, the cable 
system state vector contains a displacement value and a force value. If assumption 2 in 
section 2.1.1 were not valid and cable bending stiffness could not be neglected, the state 
vector would require four elements: displacement, slope, shear force and bending moment. 
The use of Fourier transforms 1s useful since the types of motion studied here will limited 
to free and forced harmonic vibrations. 

Within a mechanical system, the state vector at any location can be obtained by 
multiplying the state vector at some other location in the system by the transfer matrix of 
the component(s) between the two locations. In a marine cable system, the transfer 


matrices are 2 by 2 matrices represented as follows: 





a 





The elements of the transfer matrices are derived using the equations of motion. Since 
they relate the Fourier transforms of displacement and force, the transfer matrix elements 
are functions of circular frequency, @. The transfer matrix of a component ts general in 
that it does not change when the component ts installed in other systems. Cable system 
components modeled here using transfer matrices include cable segments, lumped and 
distributed masses and vibration absorbers. 


Figure 2-1 illustrates how transfer matrices are used to relate the state vectors at 


individual locations in the cable system: 


absorber 





[22] = [T1] [21] (2-3) 
[2s] = [11] [12] (21) (2-4) 
[23] = [12] [22] (2-5) 
[2(x)] = [TiCs)] [Zi] 5 x < x2 (2-6) 


Figure 2-1: Transfer Matrix Representation of Cable / Absorber System 
The transfer matrix for the uniform cable section, [T;], comes directly from the 


solution to the linear string wave equation and Is given by: 


cos(k/) 1/kT sin(k/) 
eee ior sin(k/) cos(k/) 





where: k = wave number =@/c (m’’) 

c = wave propagation speed = (T/)'* (ms) 

T =cable tension (N) 

o = cable linear mass density (kg/m) 

/ = cable length from x, to x2 (m) 
For locations on the cable to the left of x2, substitution of x for /in equation 2-7 will yield 
the transfer matrix [T,(x)] as in equation 2-6. 

The transfer matrix for the absorber, [T2], 1s given in Appendix B and tn the thesis 

by Li [5]. Appendix B also contains the transfer matrices for various vibration absorbers 
and discrete masses which may be located within a cable system. 


2.1.3 The Overall Transfer Matrix 


The product of the three transfer matrices, [T,] [T2] [Ts], equals the overall 
transfer matrix of the system. For the system shown in figure 2-1, the boundaries are 
walls. For a wall, cable displacement would be zero and cable force would be unknown. 


Thus, the state vectors at the boundaries can be represented by: 


“Live ‘ 
A= Ty’(x) (2s ) 


The overall system transfer matrix can be used to relate the two boundary state vectors: 
[Toversn] = [Ti] [T2] [Ts] (2-9) 


[24] = [Toverau] [Z1] (2219) 


_ Oe 





: I-{e <p aT 
Ty(xs) | | tar te2 | Ty’(x)) a) 


The zero displacement at the boundaries results tn: 

(ti x 0) +(tie x y (x1)) =0 (2-12) 

tl (2-13) 
This is the characteristic equation for the system. Since transfer matrix elements are 
functions of frequency, the roots of ti2 of [Toveran] are the natural frequencies of the 
system. 


The mode shapes of the system can then be plotted by setting a unity force at the wall: 


l= || Q-19 


It follows that: 
[z(x)] = [Ti0)] [zi] 5 x S x2 Colle) 
= [Ti] (T2] (T3)] (21) 5 x 2 xs (2-16) 
Note that [T3(x)] 1s expressed by substituting the quantity (x - x; ) for / into equation 2-7. 
This coordinate transformation is required because transfer matrix expressions normally 
reference a coordinate system which ts located on the component that the transfer matrix 
represents. 


In general, 


tri(x) trex) | 0 
[2(3)] = (T(9)] fa -| —— iH (2-17) 


(21) eee) 
where [T(x)] 1s the transfer matrix product of the system components between the left wall 


and x. 





The first element of [z(x)] 1s therefore: 
¥(x) = ti2(x) (2-18) 
Equation 2-18 specifies the displacement mode shape of the system. Note that t)2(x) in 
equation 2-18 is an element of the transfer matrix between the left wall and an arbitrary 
point, x, and is not the same as t)2 Of Toverat: IN equation 2-13. 
In chapter 3, this approach will be applied to a specific marine cable / absorber 


system. 





Chapter 3 


Vibration Absorber Analysis 


As mentioned in the introduction, the vibration absorber which will be the primary 


focus of this analysis is that shown below in figure 3-1: 





h=c¢+d l=e+f 


Figure 3-1: In Line Vibration Absorber 


The transfer matrix expression for this absorber was derived by Li [5] and is quite lengthy, 
requiring a page and a half to state (see Appendix B). This is because the transfer matrix 
must account for both the translational and rotational motion of the masses, damping and 
stiffness. MATLAB was used to compute the values of this transfer matrix and was used 
for all other modeling and calculation in this thesis. Appendix A contains the MATLAB 
programs STRINGI and STRING2 which calculate the transfer matrices of the absorber 


and cable segments, respectively. 
3.1 Modeling the Vibration Absorber in a Cable System 


Since the nature of this analysis was numerical, a baseline absorber system of the 
type shown in figure 3-1 was modeled. System parameters such as cable length, cable 


diameter and absorber mass were set to values which could realistically be expected in a 


17 





marine system. Optimum absorber performance can be obtained by systematically varying 


these parameters. Baseline system parameters are as follows: 


Total system length: 50.8 m 
Cable tension, T: 5000 N 
Cable density, p: | kg/m 
Absorber link mass, M: 9 kg 


Absorber link moment of inertia, I: 0.12 kg m? 


a and b dimensions: 0.03 m 

c, d, e, f dimensions: 0.2m 

Total absorber length, Lat: 0.8 m 
Stiffness, K: 10 N/m 
Damping, R: varies N/(nv/s) 


Note that the stiffness, K, is relatively insignificant. Cable tension is the dominant factor in 
providing a restoring force for the absorber. 
Absorber analysis consisted of the following steps: 
1) Obtain natural frequencies 
2) Determine damping ratio 
3) Determine displacement mode shapes 
4) Vary system parameters to optimize damping ratio 


This process 1s detailed in the following sections. 





3.2 The Cable / Absorber System With no Damping 


It is valuable to first analyze the system and demonstrate the analysis process with 
the damping value of the absorber set to zero. This keeps all calculations in the real 
domain and provides for a better understanding of the introduction of damping. 

3.2.1 System Natural Frequencies 

Chapter 2 detailed the process for determining the characteristic equation and 
natural frequencies using transfer matrices. The overall transfer matrix of the system, 
[Tovera], 1S Computed as follows: 

[Toveran] = [11] [T2] [Ts] (eo) 
where: 
[T,] and [T;] are the transfer matrices of the left and right cable segments 
as in equation 2-7 
[Tz] 1s the transfer matrix of the absorber as in Appendix B 
Due to the fixed wall boundary conditions, the system characteristic equation is 
ti9=0 (Ge 
where t)2 Is the first row, second column element of Toveran. 

A MATLAB script file, STRING4, was used to generate figure 3-2 on the 

following page, a plot of ty2 versus frequency, @. The zeroes of t)2 are the natural 


frequencies of the system. 
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T(1,2) CHARACTERISTIC FOR SYSTEM WITH ZERO DAMPING 
0.012 : 


0.01 


0.008 


0 5 10 nae? 20 25 
Figure 3-2: ty, Characteristic of Undamped System 
3.2.2 System Mode Shapes 
Once the system natural frequencies are known, transfer matrix multiplication can 
be used to determine the displacement and force mode shapes. This 1s accomplished by a 
MATLAB program called STRING6 which takes the system frequency as an input. The 


state vector at the left wall is set as follows 


0 
Zieft = 


where the first element indicates the zero displacement at the wall and the second element 


1s Set at one. This results in a mode shape which 1s normalized to a unity transverse force 


Be 





at the system boundaries. The displacement mode shape 1s calculated along the length of 
the system as described in Chapter 2: 


from left wall to absorber: 


ti ti2 | O y(x) Pe 
rotetrotea[°* Pf] es 


y(x) = ti2(x), (3-5) 
where t}2(x) 1s an element of [T)(x)] 
from absorber to right wall: 
y(x) = ti2(x), (3-6) 
where t)2(x) is an element of the matrix product [T,] [T2] [Ts(x)]. 
Similarly, the transverse force mode shape would equal the t22 element of the transfer 
matrix as a function of distance from the left wall. The resulting displacement mode 


shapes for the first three vibration modes are shown in figure 3-3 on the following page. 
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Figure 3-3; Undamped System Displacement Mode Shapes 
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3.3 Introduction of Absorber Damping 


The introduction of absorber damping takes the transfer matrix calculations into 
the complex domain. In addition to natural frequency and mode shape, damping ratio can 
now be determined. Damping ratio 1s the parameter used in determining absorber 
effectiveness. Damping ratio is used to compute the reduced damping in figure 1-2 and to 
determine if sufficient damping exists to reduce vortex-induced vibrations. 

3.3.1 System Natural Frequencies 

The same approach which was used to compute undamped natural frequencies in 
section 3.2.1 is used here to compute damped natural frequencies. As shown in section 
3.2.1, the characteristic equation of the system 1s: 

t12=0 (3-7) 
where tj2 is an element of the overall system transfer matrix. Figure 3-4 on the following 
page shows the magnitude of t;2 as a function of frequency, @. Note that when damping 
IS present, tj. is a Complex quantity. The figure also shows t)2 for the undamped system 
In the same frequency range 

The damped curve ts quite different in appearance from the undamped case. For 
alternating modes, the damped value of t;2 fails to reach a zero for vea/ values of natural 
frequencies. This behavior is significant and will later be shown to have a very physical 


explanation. 
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Figure 3-4: Magnitude(t,.) Versus Frequency for Damped vs. Undamped System 


In order to obtain all the zeroes of tj2, complex values of frequency must be considered. 
This requires the use of a computer complex root solver, MULLER (Appendix A), 
developed by Conte and de Boor [14] and converted into MATLAB format by Wilson and 
Roberts. Before proceeding further with the results of the complex root search, the 
meaning of a complex natural frequency must first be discussed. 
The harmonic motion at a point in the cable system can be described by: 
chee. (3-8) 
where A 1s displacement amplitude. 


For imaginary natural frequencies: 


© = Oreal T lela (3-9) 
Substituting: 
y(Xx) = Aekoreal + 11Mag)t (3-10) 
=A e'(oreal)t e(ounag)t (3-1 1) 


= A[cos(@,eait) + 181N(@yeait)] ee sD) 
The real part of which 1s: 
= A COS(Oreail) e tomas} (3-13) 


So, we see that for free vibrations the damped system vibrates at @,ea with a decay 


-(wimag)t 


envelope e _ The ratio of the real and imaginary parts of the natural frequency is the 
modal damping ratio: 
C a Ginger Oreal (3- 14) 


Now that the significance of a complex natural frequency 1s understood, we can 


return to the analysis of the cable / absorber system. The complex root solver, MULLER, 





is iterative and requires a starting guess for the complex frequency root. The undamped 
natural frequencies were used as initial guesses when running MULLER. With the 
assistance of Dr. Rama Rao, MIT, a routine called MULRUN (Appendix A) was 
developed which successively runs MULLER for many values of damping constant to 
obtain an output which demonstrates the optimum value of damping. Figure 3-5 below 


shows the output of MULRUN for the first vibration mode of the baseline system. 
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Figure 3-5: Demonstration of Optimum Damping Value 





The MULRUN output shows that as damping constant, R, 1s varied, there is an 
intermediate damping value which optimizes (maximizes) damping ratio. It was theorized 
that the reason damping ratio decreases at high values of damping constant is that relative 
motion between the absorber masses decreases above the optimum value. This was 
confirmed by creating a transfer matrix model of a uniform bar (see Appendix B) in place 


of the vibration absorber in the baseline system as shown in figure 3-6: 


damper 


cable 
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Uniform bar 
cable 
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Figure 3-6: Representing High Stiffness Absorber ay Uniform Bar 


The uniform bar had the same dimensions and mass as those of the absorber. It was 
shown that the cable system with the uniform bar at its center had the same natural 
frequency and mode shape as the absorber with very high damping constant. So at very 
high damping, relative rotational motion between the absorber links is restricted and the 


absorber is unable to dissipate energy. The absorber acts as a uniform bar. 
3.3.2 System Mode Shapes With Damping 


Once the system natural frequencies are known, transfer matrix multiplication can 


be used to determine the displacement and force mode shapes using the same approach as 








for the undamped case. This is accomplished by a MATLAB program called STRING6 
which takes system frequency as an Input. For the damped case, complex natural 


frequencies are input into STRING6. The state vector at the left wall 1s again: 


0 
fzal=|") (3-15) 


The displacement mode shape ts calculated along the length of the system as tn section 


3.2.2. The resulting mode shapes for the first three vibration modes are shown tn figure 


3-7, 


first mode shape 
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Figure 3-7: Damped System Displacement Mode Shapes 
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Inspection of these mode shapes reveals that for the even mode, the absorber is located at 
a modal antinode. There is no relative motion between the absorber links. The absorber 
simply rotates about its center as if it were a uniform bar. This is confirmed further by 
using MULRUN for the same modes to determine optimum damping ratio. In figure 3-8 
below, it is seen that for the even vibration mode, damping ratio is insignificant even 
though damping is present. Relative motion between the absorber links is required if the 


absorber is to dissipate energy. 
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Figure 3-8: Damping Effect in First Three Vibration Modes 





So, from initial analysis the following observations are made: 
1) An optimum value of damping constant 1s one which adds sufficient 
damping while not overly restricting motion between the two absorber 
links. 
2) Absorber effectiveness is very dependent on absorber location. 
Absorbers located at modal antinodes provide the most energy dissipation. 
These observations are consistent with the behavior of more conventional types of 
vibration absorbers, such as those installed in discrete (vice continuous) mechanical 
systems. 
In chapter 4, parameters such as absorber mass, absorber length and number of 


absorbers will be varied to study their effect on absorber performance. 
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Chapter 4 


Parameters Affecting Absorber Performance 


In previous chapters, it has been demonstrated how natural frequencies, mode 
shapes and damping ratios can be obtained for a vibration absorber in a marine cable 
system. It was shown how the damping constant, R, can be varied to arrive at an 
optimum damping ratio. The tools are now in place to study the effects of other 
parametric variations on absorber performance. Improvements in absorber performance 
will be sought by varying the following: 

1) Absorber length 
2) Absorber mass 
3) Number of absorbers installed in a single cable system 


4) Absorber location 
4.1 Required Damping 


Before proceeding with further attempts to optimize damping, it would be helpful 
to know what value of damping 1s required in the system to reduce vortex induced 
vibrations. Recall from figure 1-2 that there 1s a value of reduced damping below which 
structural damping provides no assistance 1n reducing vortex induced vibrations. 
Structures with reduced damping in this plateau range will have displacement amplitudes 
which are restricted by limit-cycle reductions to the lift forces. Thus, there is no gain in 
adding small amounts of damping which will not raise reduced damping above that value 
required to reduce the response to less than that resulting from limit cycle behavior. From 


figure 1-2, this threshold value of reduced damping 1s as follows: 








where: 

Spa ae Sein eribe (4-1) 
For the baseline system in water (see section 3.1), this results in a threshold (required) 
value of damping ratio: 

C> .04 (4-2) 
The highest value of damping ratio achieved thus far (see figure 3-8) is approximately 
C= .0012. So, we see that if the absorber being studied 1s to be of use in reducing vortex 


induced vibrations in water, further optimization is required. 


4.2 Variation of System Parameters 


4.2.1 Absorber Length 


In figure 3-7, mode shapes for the first three modes were shown. These low 
modes have very long wavelengths, resulting 1n the situation shown in figure 4-1(a) on the 
following page. When absorber length is small compared to wavelength, relative angular 
motion between the absorber links is very small, even in odd numbered vibration modes. 
Since energy dissipation relies on this relative motion, short absorber links result in poor 
absorber performance. As absorber length is increased, relative motion between the 
absorber links 1s increased, and energy dissipation increases as shown in figure 4-1(b). A 
key dimensionless parameter of use here would be the ratio of absorber length to the half 
wave length ofa given mode. This parameter is designated as the length ratio: 


LR =) Liteorber ! d (4-3) 
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Figure 4-1: Effect of Absorber Length on Damper Motion 


Variation of absorber length was accomplished in an indirect way. Instead of varying the 
absorber link length, the effect of absorber length was studied by varying A, or simply by 
comparing the damping ratios of different modes of the baseline system. Starting with 
mode 1, optimum damping values were determined for each mode. As discussed in 
section 3.3.2, even numbered vibration modes exhibited no damping effects. The best 
absorber performance (as determined by maximum damping ratio) was achieved in the 
29th mode. 
In this mode: 

hos 310) ag 

teen 

LR = .457 
In modes beyond mode 29, higher length ratio resulted in more restricted absorber motion 


and reduced damping ratios. 
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Figure 4-2 shows the damping value optimization and mode shape for mode 29. 
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Figure 4-2: Damping Optimization and Mode Shape for Mode 29 
It can be seen that with an optimum value of damping constant, damping ratio is 
approximately 1.7%. This falls well short of the value of 4% required to reduce vortex 
induced vibration amplitude in water. 


4.2.2 Absorber Mass 


In further efforts to increase the amount of damping which the absorber adds to 
the system, the mass of the absorber links was varied. A useful dimensionless parameter 
for this analysis was the mass ratio: 

MR = Mabsorber / (Latsorber Peable ) (4-4) 
This 1s the ratio of the absorber mass to the mass of the cable segment displaced by the 


absorber. 








The variation of mass ratio was performed for the 29th mode of vibration of the 
cable system. This mode was selected because it this was the mode for which optimum 
absorber length was shown. Thus the result of this variation represents an attempted 
optimization of both absorber mass and length. There is no guarantee that optimum length 
ratio remains at 0.457 when mass ratio ts varied. Figure 4-3 shows how optimum 


damping ratio varies with absorber mass: 
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Figure 4-3: Affect of Mass Ratio on Mode 29 Damping 
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Figure 4-3 shows that increasing absorber mass ratio tends to increase optimum 
damping ratio. Also, lower values of damping constant are required to achieve the 
optimum damping value when a higher mass ratio is used. This figure demonstrates that 
the baseline system (MR=22.5) has the optimum mass ratio because It results in the 
highest system damping ratio. When mass ratio was increased beyond 22.5, the system 
shifted to an even numbered vibration mode where damping is insignificant. 

If this type of analysis were to be performed on a real marine system, one would 
have to exercise care in changing absorber mass. Changing absorber mass changes the 
mass of the vibrating system and thus changes system natural frequencies. The marine 
system would be designed with some range of current velocity and corresponding vortex 
shedding frequencies in mind. If one were not careful, one could change system mass such 


that the absorber was optimized at some frequency outside the range of concern. 
4.2.3 Behavior at Realistic Values of Absorber Mass 


Admittedly, a mass ratio of MR=22.5 corresponds to a very large absorber mass, 
equaling 36% of the total cable mass. Absorbers of this extreme size were considered 
only in an attempt to achieve the required damping level of section 4.1. A more typical 
absorber mass might be no more than 10% of the cable mass. This would correspond to a 
mass ratio of MR=6.25. In this section, an absorber of this size was analyzed over many 
vibration modes and corresponding length ratios. The results are shown in figure 4-4 and 


4-5 on the following page. 
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Figure 4-4: Optimized Damping Ratio vs. Mode Number and Leneth Ratio for Mass Ratio = 6.25 
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Figure 4-5: Optimized Damping Constant vs. Mode Number and 
Leneth Ratio for Mass Ratio = 6.25 


Figures 4-4 shows that an optimum damping ratio 1s achieved in mode 53 with a length 
ratio of LR=0.83. Figure 4-5 shows how required damping decreases dramatically with 


Increasing length ratio. This behavior was discussed in section 4.2.1. 
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4.3 Number of Absorbers 


Up to this point, variation of system parameters has not resulted in achieving the 
required damping ratio of Section 4.1. It would be very convenient if one could 
multiplicatively increase damping ratio by simply increasing the number of absorbers 
installed. 

This addition of absorbers was performed for mode 3. A transfer matrix model 
was constructed for a cable system with vibration absorbers at its mode 3 antinodes. An 
absorber mass ratio of MR=6.25 was used. All other system parameters were those of the 
section 3.1 baseline system. The MATLAB programs which model this system are located 
in directory MULTABS in Appendix A. These programs are functionally identical to 
those used to evaluate the system with one absorber. 

Figure 4-4 on the following page shows the mode 3 mode shapes for the system 
with one absorber and the system with three absorbers. For both cases, a mass ratio of 
MR=6.25 and a length ratio of LR=.047 was used. The maximum damping value 


achieved in each case 1s indicated. 
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Figure 4-6: Comparison of Single Absorber to Multiple Absorbers 


Inspection of figure 4-6 reveals that optimized damping ratio roughly doubled 
when the number of absorbers was tripled. Thus, damping ratio 1s not multiplicative with 


the number of absorbers installed in the system. One should also observe from figure 4-4 








that mode 3 natural frequency decreases considerably with the additional mass loading of 


the extra absorbers. 


4.4 Absorber Location for Operation in Multiple Modes 


In previous sections, absorber performance has been looked at one mode at a time. 
Parameters such as damping constant and absorber mass have been varied and the effect 
on modal damping ratio has been observed. This has improved understanding of the 
system but it 1s not the way in which a marine vibration absorber system would be 
designed. In a system subject to vortex induced vibrations, some range of potential 
currents would exist. This range of currents would be directly correlated to an excitation 
bandwidth by the Strouhal relationship of equation 1-1. So, an absorber design problem 
would be concerned with multiple modes. This section will illustrate how the proper 
location of a single absorber can optimize performance in multiple vibration modes. 

In section 3.3.2, it was observed that the most effective location for an absorber 
was at a modal antinode. For mode 3, these locations would be L/6, L/2, and 5L/6. Now 
suppose we wanted to design an absorber system which would be effective for modes 3 
and 4. The mode 4 antinodes are at L/8, 3L/8, 5L/8. The two modes share no common 
antinodes. In fact, the L/2 antinode for mode 3 1s a mode 4 node. So, an absorber located 
at L/4 offers no absorber effect for mode 4. 

A compromise was arrived at by constructing a model of a system with an 
absorber at 7L/16. This location is half way between the 3L/8 mode 4 antinode and the 
L/2 mode 3 antinode. The absorber used had a mass ratio of MR=6.25. The MATLAB 


programs for modeling and analyzing this system are located in directory OFFCENTER in 
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Appendix A. Table 4-1 below compares absorber effectiveness for the compromise 


absorber location versus absorber locations optimized for a single mode: 


a 
Effectiveness Mode Effectiveness 

Ratio 

- nll 
Ratio 


Table 4-1; Dual Mode Absorber Optimization 







The results indicate that relocation of the absorber causes some loss of mode 3 
damping ratio while raising mode 4 damping ratio from zero to a value equal to that of 
mode 3. This is the type of tradeoff would have to be performed in the design of a real 
marine cable absorber system. 

Figure 4-7 on the following page shows the mode 3 and 4 mode shapes with the 
absorber located at 7L/16. This figure demonstrates the utility of the STRING6 mode 
plotter in visualizing unorthodox mode shapes. 

So, absorbers can be located for effectiveness in both even and odd numbered 
modes. In choosing absorber location, the designer of a marine cable system would have 
to consider the excitation bandwidth and vibration modes within that bandwidth. Ifa 
single absorber location could not provide energy dissipation tn all possible excitation 


modes, then multiple absorbers would need to be installed. 
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Figure 4-7: Mode 3 and 4 Shapes with Absorber at 7L/16 





4.5 Consideration of a Different Absorber Type 


Up to this point in the analysis, considerable insight has been gained in the 
behavior of the hinged, in-line, vibration absorber in the marine cable system. However, 
this type of absorber has not been able to produce the necessary damping ratio required to 
reduce vortex induced vibrations in a cable system in water. (See section 4.1). 

In this section, a simpler and more conventional absorber will be analyzed. This 


absorber is shown in figure 4-8 below: 


cable 





Damper, R 





Spring, k 
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Figure 4-8: Hanging-Miass-Spring Dashpot Absorber 
The absorber is of the conventional mass-spring-dashpot type. It is usually associated 
with discrete mechanical systems but could also be used to dissipate energy in a cable 


system. Li and Vandiver [10] derived the following transfer matrix for this absorber: 


l 0 
[Tl=| w*mGo@R+tk) ay, 
k - om + j@R 





A numerical model of a cable system using this absorber was constructed in 
MATLAB. The cable itselfis identical to that of the baseline system of section 3.1. 
However, the absorber of figure 4-8 is located at the center of the cable in place of the in- 
line bending type absorber. The programs comprising this model are located in directory 
HANG in Appendix A. 

Den Hartog [1] formulated a method for choosing optimum absorber mass and 
stiffness for this type of absorber. His formulation was for discrete systems but modal 
analysis of the continuous cable system allows Den Hartog’s method to be applied in an 
approximate way. The approach uses the following parameters: 


tt = m/ M = mass ratio = absorber mass/ main mass 
o = V¥k/m =natural frequency of absorber 


€)= VK /M =natural frequency of main system (cable without absorber) 

f= /Q = frequency ratio 

Den Hartog determined that the correct tuning of the mass spring dashpot 
absorber results from the frequency ratio being set as follows: 

f=1/(1 +n) (4-6) 

This equation was used to set absorber stiffness. To do this, an example mass ratio of 
[t = 0.2 was selected. The main mass was set at the modal mass of the cable which equals 
one half of total cable mass. © was set at the cable first mode natural frequency. 


Equation 4-6 then allowed computation of the absorber stiffness, k. 
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The following values were used in the analysis of the vibration absorber for the first 


vibration mode of the cable system: 
m = 5 kg = absorber mass 
M = 25 kg = cable modal mass = % of total cable mass 
u=m/M=0.2 


C) = 4.44 rad/sec = cable first mode natural frequency without absorber 


Using equation 4-6: 


I k 
@ = 3.70 rad/sec = ,/— 
m 


k = 68.5 N/m, the correct stiffness for optimum tuning 
Once stiffness and mass were set, MULRUN was used to determine the optimum 
damping constant for the absorber designed tn the procedure above. Figure 


4-9 below shows the resulting damping ratio performance and first mode shape: 
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Figure 4-9: DamperOptimization for First Mode with Hanging Mass Spring Dashpot Absorber 





It can be seen that optimum damping ratio 1s approximately 0.22, far higher than 
any damping ratio achieved with the in line absorber. This damping ratio is also above the 
critical value of € = .04 necessary to reduce vortex induced vibration amplitude in water. 
This result also provides a reasonable reconciliation of results between the Den Hartog 
method and the transfer matrix method used in this thesis. Den Hartog states that the 
optimum damping constant for the mass spring dashpot absorber can be determined by the 
following formula: 


Roptimum 
2m 


= ee are 
=[3u/80+y))] C_ 
For a mass ratio of tt = 0.2, this results in an Roptimum Of 9.25 N/(m/s). This corresponds 
roughly with the optimum damping constant of 15 N/(m/s) shown in figure 4-7. 

The hanging mass-spring-dashpot absorber represents a large improvement in 
absorber effectiveness over the bending type absorber discussed elsewhere in this thesis. 


This absorber has some possibility for use in marine cable systems and further analysis 


Should be conducted. 
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Chapter 5 


Conclusion 


In this thesis, the transfer matrix method was used to construct various numerical 


models of marine cable systems. The performance of various vibration absorber 


configurations was studied. The following are the main findings of this research: 


l) 


2) 


Introduction of damping in the transfer matrix model causes the system 
characteristic equation to have complex roots. The use of a computer complex 
root solver was introduced and was shown to be useful in calculating system 
natural frequencies and damping ratios. 

Upon understanding the significance of complex natural frequencies, computer 
programs using transfer matrices were created to plot system mode shapes. 
The ability to visualize mode shapes was very important in understanding 
system behavior. 

The use of an in-line bending type vibration absorber was studied in some 
depth. A method was devised for tuning such an absorber to achieve 
maximum damping within a marine cable system. Variation of parameters such 
as absorber mass, length and number of absorbers was performed to study the 
effect of these variations on absorber performance. Even after optimization, 
the in-line absorber was unable to demonstrate sufficient damping to reduce 


vortex-induced vibrations in water. 
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4) The use of a hanging mass-spring-dashpot absorber was also studied for use in 
a marine cable system. This absorber was shown to have the ability to create 
sufficient damping to reduce vortex induced vibrations in water. 

The main focus of this thesis was to study the performance of the in-line bending 
type absorber in a finite length marine cable system. Although it was shown that this 
absorber provides insignificant damping to reduce vortex-induced vibrations in water, 
considerable insight was gained in this thesis. A working model was created which has the 
ability to calculate system natural frequencies, mode shapes and damping ratios. 
Additionally, it should be noted that the in-line absorber could be useful in reducing vortex 
induced vibrations in air. The studies in this thesis would be useful in optimizing the 
absorber for this application. 

Areas for potential research related to this thesis include: 

1) Use of the transfer matrix models in this thesis to study other marine cable 
system characteristics, such as the effects of distributed masses within the 
system. Li [5] performed considerable work in this area but the ability to 
visualize mode shapes would provide a good supplement his work. 

2) Extension of the analysis procedure in this thesis to other systems which can be 
represented by transfer matrices. 

3) Further study of the mass-spring-daspot absorber for use in marine cable 


systems. 
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Appendix A: MATLAB Files 


This appendix contains all MATLAB files used in this thesis. The files are divided 


into four directories, each corresponding to a different cable system model. The 


following is a list of these directories and a description of the system they model: 


CENTER: A cable system 50 meters long with an in line bending type absorber at 


its center. 


OFFCENTER: A cable system 50 meters long with an in line bending type 


absorber located at 7L/16. 


MULTABS: A cable system 50 meters long with three in line bending type 


absorbers, each located at a mode 3 antinode. 


HANG: A cable system 50 meters long with a hanging mass spring dashpot 


absorber located at its center. 


Variables and Dimensions 


The following variables are used throughout the MATLAB programs in this 


Appendix: 


Variable 
a, b, ema. ec, f 
M1, M2 


Corresponding Parameter 


cable dimensions 

absorber link masses 

cable tension 

stiffness 

damping constant 

absorber link moments of inertia 
cable mass density 

cable wave propagation speed 
cable wave number 

circular frequency 


50 








CENTER 


The following 1s a typical sequence for using the model CENTER to determine 


system natural frequencies and mode shapes and optimize system damping: 


1) In MATLAB, change directory to CENTER. 
2) Edit STRINGI for desired absorber properties and STRING2 for desired cable 
properties. 
3) In STRING4 and STRINGI, ensure damping constant Is set to zero. Run STRING4 
to plot t12 of [Toveran] versus frequency. The zeroes of this plot are the undamped system 
natural frequencies. 
4) If desired, plot undamped system mode shapes as follows: 
a) For mode of interest, enter natural frequency in MATLAB window. 
b) Run STRING6 to plot mode shape. 
c) Repeat for other modes of interest. 
5) Perform absorber optimization for individual modes as follows: 
a) Edit MULRUN to iterate over desired range of damping constants. 
b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 
6) Set damping constant to global in STRINGI. Run MULRUN to plot real and 
Imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 
7) Plot damped system mode shape as follows: 
a) Set damping constant to the optimum value in STRINGI. 
b) Enter the complex natural frequency in the MATLAB window. 
c) Run STRING6 to plot mode shape. 


8) Repeat steps 5 through 7 for other modes of interest. 


2) 





%o STRING1 


% This file contains the transfer matrix for a vibration 
% absorber. 


% Initial data 


T=5000; 
a= 03: 
b=03- 
c=: 
G=.2: 
c= 
f=2' 
M1=9.0; 
M2=M 1: 
Hei: 
Oa ee 
Vall: 
K= 10: 
=4000; 
h=c+d; 
l=e+f; 


al=(T*d)+(H*c)+(K*a*2)-(w*2*11)+(i*w*R*b’2); 


a2=(T*f)+(H*e)+(K*a°2)-(w°2*12)+(1* w*R*b°2); 
a3=(K*a 2)+0*¥w*R*b 2); 


g1=(-(M1*1+M2*f)*al- 
(M1 *c*a3)+(w°2*M 1 *M2*£*c°2))/((M2 *e*al)+(M 1 *d+M2*h)*a3+(w°2*M1*M2*c*e*d)); 


g2=(-(M2*f*h*2)- 
(M1*1*d°2)+(a1*1+a3*h)/w°2)/((M2*e*a1)+(M 1 *d+M2*h)*a3+(w°2*M 1*M2*c*e*d)); 


a 
((M1*1*d*c)-+(al *l+a3*h)/w*2)/((M2*e*a 1 )-+(M 1*d-+M2*h)*a3+(w*2*M 1*M2*c*e*d)): 


g4=-((M1*14+M2*f)*a3+(M 1 *c*a2)+(w°2*M 1*M2*e*c*f))/((a2*h+a3*1)/w2- 
(MT=d*F2-+-M2*t 2*h)); 


g5=((a2*h+a3*1)/w°2+(M2*e*f*h))/((a2*h+a3*1)/w*2-(M 1 *d*1°2+M2* f° 2*h)); 


g6=(-a2*(M1*d+M2*h)-(M2*e*a3)-+(w72*M 1*M2*e"2*d))/((a2*h+a3*1)/w72- 
(M1*d*1°2+M2*f°2*h)): 











t(1, 1)=gl+g3*((g4+g1*g6)/( 1-g3*g6)); 
t(1,2)=g2+g3*((g5+g2*g6)/(1-g3*g6)); 
t(2,1)=(g4+g1*g6)/(1-g3*g6); 
t(2,2)=(g5+g2*g6)/(1-g3*g6); 





To STRING2 
% This file calculates the transfer matrix for a string of length L. 
%oInitial data 


= 5000; 
[=25, 
p=i0: 
c=sqrt(T/p); 
keyve: 
ere =cos(k*L). 
tele =sintk (k= 1): 
(A= -k*T*sin(k*L); 
ie —cos( k=): 





Jo STRING4 


% This file plots the overall T(1,2) for a system with 3 components. 
W=0; 
TrotlZ=0: 


Ko: 

fomn= 12250: 
w=.1*n; 

W(nj=w; 

string? 

lee 

jis 

string 1 

=e 

Tto=Tl* 127 is: 
TtotlZan)= Ten, 2): 
end 
plot(W,abs(Ttot12)) 
grid 

xlabel(’ w(rad/sec)’) 
ylabel( T(1,2)’) 








Zo STRINGS 


% This file contains the transfer matrix for a uniform bar. 
% Initial data: 


% bar length 
p="6" 


M=18.0; 
T=5000; 
IbeM*Ib*w°2/T: 


k=1/(2*(6+Ib)): 
t(1,1)=k*4*(3-Ib); 
t(1,2)=k*12*Ib/T: 
t(2,1)=-k*M*w’2*(12-Ib): 
t(2,2)=k*4*(3-Ib): 





Jo STRING6 


Jo (ES Rie Peers THE MODE SHAPES OF AS TRING 
Jo WITH A DAMPER AT ITS CENTER. BEFORE RUNNING, 
Jo ENTER FREQUENCY IN THE MATLAB WINDOW AND 
% ENTER DAMPING IN STRINGI. 

y=0; 

08 

T=5000; 

p=0, 

c=saqrt(T/p); 

kK=Wice 

xO oe). o 


jfoNn=12251 

tile =cost ke x(n): 
t1(1,2)=sin(k*x(n))(k*T): 
tl(2,1)=-k*T*sin(k*x(n)); 
tl(2,2)=cos(k*x(n)); 
y(n)=tl(1,2); 

f(n)=th(2,2); 

end 


string 1 
isiei=cOs( 2a kh) 

Gs 2) =sin(2o7k)/(K* 1). 
iG 292) =—COs| 255k). 

6. =-k* T7sin( 25k): 


for n=258:509 

tne =costie=(xim)=25.0))- 
tr(1,2)=sin(k*(x(n)-25.8))/(k*T); 
the 2e =" Sime (x (m)- 25.8): 
Uiieee = COS( xia ora): 
Seige 

y(n)=tf(1,2); 

f= 22) 

end 


plot(x(1:251),real(y(1:251)),x(258:509),real(y(258:509))) 


grid 
xlabel(’x location’) 


a 





ylabel(’real(displacement mode shape)’) 
title? 3rd MODE MODE SHAPE’) 
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Jo STRINGI4 
% This file computes T(1,2) for a system with 3 components. 


function T12=string14(w) 
string2 

PSE 

T3=t; 

string 1 

t2=t: 

Ttot=T1*T2*T3; 
T12=Ttot(1,2); 





% MULLER 


function [zros]=muller(fn,nprev,np,maxit,ep1,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,ep1,ep2,fnreal,zeros) 

% This function is called by the driver program mulrun.m 

% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, 1.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 

% by S. D. Conte and Carl de Boor, "Elementary Numerical 

% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

% 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
% no better estimates are available. 

% ep|: relative error tolerance on the root. 

% ep2: error tolerance on the function. 

Y% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% fnreal: a logical number; if TRUE, the program forces all approximations to 


% all the roots to be real. This makes it possible to use this routine 

J even if f(x) is defined only for real x. 

% Here we can take fnreal=1, if we only want to fine the real roots; we 
% can take fnreal=0 if we want to find both the real and complex roots. 


fn=’string14’; 
epsl=max(epI,1.e-12); eps2=max(ep2,1.e-20); 
for i=snprev+l:np; kount=0; 

while | 


% Compute the first three estimates for zero as 
% zero, zerot+0.5, zero-0.5 


zero=zros(1); h=0.00001+j*0.00001; zp=zero+h; zm=zero-h; 
[fzr,dvdfip,kount,zros]=dflate(fn,zp,i,kount,zros); 
[fzr,fzrprv,kount,zros]=dflate(fn,zm,1,kount,zros); 
hprev=-2.0*h; dvdf1p=(fzrprv-dvdf1p)/hprev; 
(fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 


Y% Do while kount<maxit or his relatively big or fzr=f(zero) 
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ee 


% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 


divdfl=(fzrdfil-fzrprv)/h; divdf2=(divdf1-dvdf1p)/(h+hprev); 
hprev=h; dvdflp=divdf1; c=divdfl+h*divdf2; 
Sei—C ¢-4.0 “izrdi divdt2: 
if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<O, den=c-sar; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfi/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 
(fzr,fzrdfi,kount,zros ]=dflate(fn,zero,i,kount,zros); 
i zios )=—zer0 1.001, break. end 
Jo Check for convergence 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
To Check for divergence 
if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end; % end kount <= maxit 
iPZ20Os( =—ZerOeuO | break end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
ZiOS(1)=Zero: 
end; % for loop 
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%o MULRUN 


x=0; 
y=0; 
global R 
% input specs for Muller 
torn 5 
R= 106 *(n- 1); 
KG=kK: 
% tolerance, name of function that computes T12 etc: 
mae Nee: % Starting guess for the root, can be real or complex 
fn = ’string14’; % name of program that computes T12 
maxit = 500: % max # of iterations 
epl = _ le-5; %tolerance on the root (start with 
% values, le-2or so.) 
tolfcn = le-5; % relative tolerance on T12 
fnreal = 0; Ja— complex search, — 1, real seatem 
nprev = 0; % leave these as 0 
np= 1; % and | 


ep2 = tolfcn*string14(zrs); % This sets the tolerance on T12 to be 
% le-5 of 712 estimated at your 
Jo guess 


% Call Muller 


root = muller(fn,nprev,np,maxit,ep1,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 


subplot(311) 
plot(x,imag(y)) 

grid 

title’ 4th mode’) 
Xlabel(’R’) 
ylabel(’imag(nat freq)’) 
subplot(3 12) 
plot(x,real(y)) 
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lower 


Starting 





grid 

xlabel(’R’) 

ylabel(’ real(nat freq)’ ) 
subplot(3 13) 

plot(x,zeta) 

grid 

xlabel(’R (N/(m/s))’) 

ylabel(’ modal damping ratio’) 
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To DEE ATE 


% This is the subroutine needed by muller. 

function [fzero,fzrdfi,kount,zros]=dflate(F,zero,1,kount,zros); 

% (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,1,kount,zeros) 

% This function which ts called by function muller is translated 

% from the Fortran routine given by S. D. Conte and Car! de Boor, 

% “Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 


kount=kount+1; fzero=feval(F,zero); fzrdfi=fzero; 
Neleareliin: enc 
forj—2"1 
den=zero-zros(j-1); 
if abs(den)==0.0 
zros(i)=zero* 1.001; return 
else 
fzrdfl=fzrdfi/den; 
end 
end 
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OFFCENTER 
The following is a typical sequence for using the model OFFCENTER to 


determine system natural frequencies and mode shapes and optimize system damping: 


1) In MATLAB, change directory to OFFCENTER. 
2) Edit STRINGIA for desired absorber properties and STRING2AL and STRING2AR 
for desired cable properties. 
3) In STRING4A and STRINGIA, ensure damping constant Is set to zero. Run 
STRINGA4A to plot ty2 Of [Toveran] versus frequency. The zeroes of this plot are the 
undamped system natural frequencies. 
4) If desired, plot undamped system mode shapes as follows: 
a) For mode of interest, enter natural frequency in MATLAB window. 
b) Run STRING6A to plot mode shape. 
c) Repeat for other modes of interest. 
5) Perform absorber optimization for individual modes as follows: 
a) Edit MULRUN to iterate over desired range of damping constants. 
b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 
6) Set damping constant to global in STRINGIA. Run MULRUN to plot real and 
Imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 
7) Plot damped system mode shape as follows: 
a) Set damping constant to the optimum value in STRINGIA. 
b) Enter the complex natural frequency in the MATLAB window. 
c) Run STRING6GA to plot mode shape. 


8) Repeat steps 5 through 7 for other modes of interest. 





Jo STRINGIA 


% This file contains the transfer matrix for a vibration 
% absorber. 


% Initial data 


T=5000; 


K=10; 

R=100000; 

h=c+d; 

l=e+f; 
al=(T*d)+(H*c)+(K*a°2)-(w°2*11)+C*w*R*b°2); 


a2=(T*f)+(H*e)+(K*a°2)-(w*2*12)+(i*w*R*b°2); 
a3=(K*a*2)+(i*w*R*b°2); 


g1=(-(M1*1+M2*f)*al- 
(M1*c*a3)+(w°2*M 1*M2*f*c°2))/((M2*e*al )+(M 1 *d+M2*h)*a3+(w°2*M1*M2*c*e*d)); 


g2=(-(M2*f*h*2)- 
(M1 *1*d*2)+(al*1+a3*h)/w72)/((M2*e*al)+(M 1 *d+M2*h)*a3+(w72*M 1*M2*c*e*d)); 


gs=- 
((M1*1*d*c)+(al *1+a3*h)/w*2)/((M2*e*al)+(M 1 *d+M2*h)*a3+(w°2*M | *M2*c*e*d)); 


g4=-((M1*14+M2*f)*a3+(M 1 *c*a2)+(w°2*M 1*M2*e*c*f))/((a2*h+a3*])/w* 2- 
(M 1*d*1°2+M2*f2*h)); 


g5=((a2*h+a3*1)/w°2+(M2*e* f*h))/((a2*h+a3*1/w°2-(M 1 *d*1°2+M2*f°2*h)); 


g6=(-a2*(M1*d-+M2*h)-(M2*e*a3)-+(w*2*M 1*M2*e72*d))/((a2*h-+a3*1)/w72- 
(M1 *d*1°2+M2*f2*h)); 
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t(1,l)=g1+¢3*((¢4+g1*26)/(1-g3*g6)); 
t(1,2)=g2+g3*((g5+g2*g6)/(1-g3*g6)); 
t(2,1)=(g4+g1*g6)/(1-g3*g6); 
t(2,2)=(g5+g2*g6)/(1-g3*g6); 
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% STRING2ZAL 


%* This file calculates the transfer matrix for a string of 


Penegt huwk, 
% located on the left side of the absorber 


Initial data 

T=5000; 

G—(7/ toy 50 - 

p=1.0; 

S=SCre(T/ pie 

K=w/c; 

Pa, Vi =eos ties Ty): 5 

ol ysl be, (KT) > 
m2 yak ~ T* Sag k= 1) > 
2. 2h COS tk il): 
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Jo STRING2ZAR 


% This file calculates the transfer matrix for a string of length L, 
% located on the right side of the absorber. 


% Initial data 


P=5000: 

L=(9/16)*50; 

p=1.0; 

c=sqrt(T/p); 

K=Ww/c: 

ee =cos(k°L): 

tem j=sin(k Lj k* 1): 
eka Sint k 12): 
ties) =COs(K" LL); 
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%o STRING4A 


% This file plots the overall T(1,2) for a system with 3 components. 
=O. 
Ttot12=0; 


i. 

fer n=1:200; 

=, lie 

W(n)=w; 

string2al 

Wiles 

string2ar 

TSS 

stringla 

Se 

Tioetia ters: 
MotlZ(n)= tow ),2): 
end 
plot(W,abs(Ttot12)) 
grid 

xlabel(’ w(rad/sec)’) 
ylabel(’T(1,2)’) 
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Jo STRING6A 


Jo Tite EeOtS Min viG De Snares Or A STRING 
Jo WITH AN OFF CENTER DAMPER. ENTER FREQ. 
%o FROM THE MATLAB WINDOW AND DAMPING 
To IN STRINGIA 

x=0; 

ye 

T=5000; 

p=1.0; 

c=sqrt(T/p); 

k=wse: 

x=0:.1:50.8; 


fonn= 12220 

BiG =costk extn) ): 

tl(] ,2)=sin(k*x(n))/(k*T); 
tl(2,1)=-k*T*sin(k*x(n)); 
tl(2,2)=cos(k*x(n)); 
y(n)=tl(1,2); 

hn) =t( 2.2); 

end 


Stringla 
ts(1,1)=cos(21.9*k); 
SG) =sin(Z21.9*k (k= 1). 
es(2. 2) =Ccos(2 197k): 
is(Z=-k* 1 since 97k), 


for n=227:509 

tml cost (xti)-22.7)): 

tiles J=SIn(kexieee, 2 yy (ke Ll): 
the =-k PP sin( ke 0x m)-22.7)): 
th(2, 2 j—COs(k (x ()=22-7 )); 

iets tty: 

y(n)=tf(1,2); 

hn) =the2 2): 

end 


plot(x(1:220),real(y(1:220)),x(227:509),real(y(227:509))) 


grid 
xlabel(’x location’) 
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‘ 
ho 
‘ 


ylabel(’real(displacement mode shape)’) 
title’’4TH MODE MODE SHAPE, DAMPER AT 3L/8’) 





%o STRINGI4A 


% This file computes T(1,2) for a cable with an off center absorber. 


function T12=string14a(w) 
Sstring2al 

Bist 

Sstring2ar 

T3=t; 

Stringla 

T2=t; 

Ttot=T1*92*T3; 
tion 2), 





To MULLER 


function [zros]=muller(fn,nprev,np,maxit,ep! ,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,ep!,ep2,fnreal,zeros) 

% This function is called by the driver program mulrun.m 

% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, 1.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 

% by S. D. Conte and Carl de Boor, “Elementary Numerical 

% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
%o Mechanics Department, University of Alabama. 

Jo 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
Jo no better estimates are available. 

% ep|: relative error tolerance on the root. 

% ep2: error tolerance on the function. 

% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% f{nreal: a logical number; if TRUE, the program forces all approximations to 


Jo all the roots to be real. This makes it possible to use this routine 

Yo even if f(x) 1s defined only for real x. 

Zo Here we can take fnreal=1, 1f we only want to fine the real roots; we 
% can take fnreal=O if we want to find both the real and complex roots. 


fn=’string 14a’; 
eps l=max(ep1,1.e-12); eps2=max(ep2, 1 .e-20); 
for i=nprev+1:np; kount=0; 

while 1 


% Compute the first three estimates for zero as 
% zero, zero+0.5, zero-0.5 


zero=zros(1); h=0.0000 1+j*0.00001; zp=zero+h; zm=zero-h; 
[fzr,dvdf1p,kount,zros]=dflate(fn,zp,i,kount,zros); 
[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 
hprev=-2.0*h; dvdfl p=(fzrprv-dvdf1p)/hprev; 
[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 


% Do while kount<maxit orh is relatively big or fzr=f(zero) 
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Ru: 


% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 


divdfl=(fzrdfi-fzrprv)/h; divdf2=(divdfl-dvdf1p)/(h+hprev); 
hprev=h; dvdflp=divdfl; c=divdf1l+h*divdf2; 
SOr—cec siz “divalz: 
if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<O, den=c-saqr; 
Clce,dem—c tsar. end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfi/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 
(fzr,fzrdfi,kount,zros]=dflate(fn,zero,i,kount,zros); 
if ze@s(i)==zero* | OOi break: end 
% Check for convergence 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
Jo Check for divergence 
if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end co Cha skount <— maxit 
if 2test)j==zcr0- 1 O01, break end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
POS) —Zere, 
end; % for loop 


a3 








Jo MULRUN 


Of 2H HR SR ARR OR SR A AR RCC A A OO KK A HK 


x=0; 
y=0; 
global R 
% input specs for Muller 
for n=1;21 
R=10000*(n-1); 
ICO EIS 
% tolerance, name of function that computes T12 etc: 
ES OS. % Starting guess for the root, can be real or complex 
fn =’stringl4a°; | % name of program that computes T12 
Mmaxit = 500. % max # of iterations 
epl = le-5; Ytolerance on the’ root’ (start with 
% values, le-2or so.) 
tolfcn = le-5; % relative tolerance on T12 
itiecieaU, % = 0 complex search, = 1, real search 
prey = 0; % leave these as 0 
np = 1; Jo and | 
ep2 = tolfcn*string 14a(zrs); % This sets the tolerance on T12 to be 
YJ \le-> of T12~ estimated at your 
Jo guess 


% Call Muller 


root = muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 


subplot(311) 
plot(x,imag(y)) 

grid 

title’ 4th mode’) 
xlabel(’R’) 
ylabel(’imag(nat freq)’) 
subplot(312) 
plot(x,real(y)) 
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a 





grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 
subplot(3 13) 

plot(x,zeta) 

grid 

xlabel(’R’) 

ylabel(’ modal damping ratio’) 


7) 


— 





Jo DFLATE 


% This is the subroutine needed by muller. 

function [fzero,fzrdfl,kount,zros ]=dflate(F,zero,1,kount,zros); 

% (*fzero,fzrdfl,kount,zeros* )=dflate(F,zero,i,kount,zeros) 

% This function which is called by function muller is translated 

% from the Fortran routine given by S. D. Conte and Carl de Boor, 

% "Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 


kount=kount+1; fzero=feval(F,zero); fzrdfl=fzero; 
if i<2, return; end 
for J=2:1 
den=zero-zros(j- 1); 
if abs(den)==0.0 
@esij=zero |.OU1- return 
else 
fzrdfl=fzrdfl/den; 
end 
end 
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MULTABS 


The following is a typical sequence for using the model MULTABS to determine 


system natural frequencies and mode shapes and optimize system damping: 


1) In MATLAB, change directory to MULTABS. 
2) Edit STRINGIA for desired absorber properties and STRING2AI and STRING2AO 
for desired cable properties. 
3) In STRING4A and STRINGIA, ensure damping constant Is set to zero. Run 
STRING4A to plot ti2 of [Toveran] versus frequency. The zeroes of this plot are the 
undamped system natural frequencies. 
4) If desired, plot undamped system mode shapes as follows: 
a) For mode of interest, enter natural frequency in MATLAB window. 
b) Run STRING6A to plot mode shape. 
c) Repeat for other modes of interest. 
5) Perform absorber optimization for individual modes as follows: 
a) Edit MULRUN to tterate over desired range of damping constants. 
b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 
6) Set damping constant to global in STRINGIA. Run MULRUN to plot real and 
imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 
7) Plot damped system mode shape as follows: 
a) Set damping constant to the optimum value in STRINGIA. 
b) Enter the complex natural frequency in the MATLAB window. 
c) Run STRING6A to plot mode shape. 


8) Repeat steps 5 through 7 for other modes of interest. 





%o STRINGIA 


% This file contains the transfer matrix for a vibration 
% absorber. 


% Initial data 


T=500G: 
=O) 
b=.03: 
c=. 
d= )- 
e=.2: 

f= 2: 
M1=9; 
M2=M1; 
H=5000; 
leo 2: 
Pere 
K=10; 
global R; 
h=c+d; 
l=et+f; 


al=(T*d)+(H*c)+(K*a°2)-(w°2*11)+(C*w*R*b°2); 


a2=(T *f)+(H*e)+(K*a°2)-(w*2*12)+(i*w*R*b*2); 
a3=(K*a°2)+(i*w*R*b°2); 


g1=(-(M1*1+M2*f)*al- 
(M1*c*a3)+(w*2*M1*M2*f*c72))/((M2*e*al )+(M 1 *d+M2*h)*a3+(w*2*M1*M2*c*e*d)); 


g2=(-(M2*f*h72)- 
(M1 *1*d72)+(al *l-+4a3*h)/w72)/((M2*e*al)-+(M 1 *d+M2*h)*a3+(w72*M1*M2*c*e*d)); 


gi=- 
((M1*1*d*c)+(al *1+a3*h)/w"2)/((M2*e*al)+(M 1 *d+M2*h)*a3+(w*2*M 1*M2*c*e*d)); 


g4=-((M1*1+M2*f)*a3+(M1*c*a2)-+(w*2*M 1*M2*e*c*f))/((a2*h-+a3*l)/w*2- 
(M1*d*I°2+M2*f-2*h)); 


g5=((a2*h+a3*l)/w*2+(M2*e*f*h))/((a2*h+a3*1)/w2-(M1*d*1°2+M2*f-2*h)); 


g6=(-a2*(M1*d+M2*h)-(M2*e*a3)+(w°2*M 1*M2*e°2*d))/((a2*h+a3*1)/w2- 
(M 1*d*1°2+M2*f2*h)); 


SO 








t(1,lj=gl+g3*((g4+g1*g6)/(1-g3*g6)); 
t(1,2)=g2+g3*((g5+g2*g6)/(1-g3*g6)); 
t(2,1)=(g4+g1*g6)/(1-g3*g6); 
t(2,2)=(g5+g2*g6)/(1-g3*g6); 


S| 








% STRING2AI 


% This file calculates the transfer matrix for a string of length L/3, 
% located on the left side of the abosorber 


% Initial data 

T=5000; 

L=(148)*50; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

Tel) —Cos(k le): 

(Ce? j=saatk*Ly(k*T): 
ee) sintk LL): 
He 2 =cos(ikew): 








% STRING2AO 


% This file calculates the transfer matrix for a string of length L/6, 
% located on the left side of the abosorber 


% Initial data 

T=5000: 

P=C GO) 50; 

p=1.0; 

c=sqrt(T/p); 

kw: 
cGap=cos(k*L): 
(enim hike TD): 
U2 kT *sin(k*L); 
2 )—COstk IL); 





%o STRING4A 


% This file plots the overall T(1,2) for a system with absorbers located 
% at the mode 3 antinodes. Before running, set R to the desired value 
% in STRINGIA. 


W=0; 
Ter 1 2=0; 


Ru; 

foen=1- 100: 
w=100+.1*n; 
W(n)=w; 

string2ao 

Tht 

string2al 

Yee: 

Stringla 

Pet: 

Memeo ls-i2*T3*T2"T |; 
Trot 1 2(m)=T tot( 1,2); 
end 
plot(W,abs(Ttot12)) 
grid | 

xlabel(’ w(rad/sec)’) 
ylabel(’T(1,2)’) 
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op STRING6A 


% THISMICE PLOTS THE MODE SHAPES OF A STRING 
Jo WITH ABSORBERS AT THE MODE 3 ANTINODES. 
% Nek WiEINC y 

Jo FROM THE MATLAB WINDOW AND DAMPING 

Jo IN STRINGIA 


p=1.0; 
c=sqrt(T/p); 
Kaw ce; 
x=0:.1:52.4; 


for n=1:84 

tlhGi I =costk* x(n): 

tl(1 ,2)=sin(k*x(n))/(k*T); 
th(2, 1)=-k*T*sin(k*x(n)); 
tl(2,2)=cos(k*x(n)); 
y(n)=tl(1,2); 

f(n)=tl(2,2); 

end 


stringla 


ts@@l, [J=cos(3.33 98); 
(GCI) —sinte 23° (k@)): 
tso(2,2)=cos(8.33*k); 
tso(2, 1 )=-k*T*sin(8.33*k); 


tsi(1,1)=cos(16.67*k); 
tsi(1,2)=sin(16.67*k)/(k*T); 
tsit2-2 )=cas(l G67 hk): 
tsi(Z, | =-k* [*sin(16.67*k); 


for n=92:259 

tle(1, 1)=cos(k*(x(n)-9.1)); 
tle(1,2)=sin(k*(x(n)-9. 1))/(k*T); 
the(2,1)=-k*T*sin(k*(x(n)-9.1)); 
tle(2,2)=cos(k*(x(n)-9.1)); 
(co fe. 

y(n)=tf(1,2): 





iim) =th( 2:2); 
end 


fon m=207 434 
tre(1,1)=cos(k*(x(n)-26.6)); 
tre(1,2)=sin(k*(x(n)-26.6))/(k*T); 
tre(2, 1 )=-k*T *sin(k *(x(n)-26.6)); 
trc(2,2)=cos(k*(x(n)-26.6)); 
tf=tso*t*tsi*t*tre; 

y(n)=tf(1,2); 

ft 2,2): 

end 


for n=442:525 

tr(1,1)=cos(k*(x(n)-44.1)); 

2) =sin( kn) 44 (kn); 

tr(2,1)=-k*T*sin(k*(x(n)-44. 1)); 

tr(2,2)=cos(k*(x(n)-44.1)); 

HetsOol sit tcl t th, 

y(n)=tf(1,2); 

fm=ti2 2) 

end 

plot(x(1:84),real(y(1:84)),x(92:259), real(y(92:259)),x(267:434),real(y(267:434)),x(442:525), real(y(442:5 
grid 

xlabel(’x location’) 

ylabel( real(displacement mode shape)’) 

title’ 3rd MODE MODE SHAPE, DAMPERS AT ANTINODES’) 
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%o STRINGI4A 


% This file computes T(1,2) for a cable with absorbers located at the 
% mode 3 antinodes. 


function T12=string14a(w) 
string2ao 

tet 

String2ai 

T=t 

Stringla 

ae =t: 


MméoHT 1 *heei3*T2*13* Tae; 
SU ZeihieiGl oe 
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Jo MITE EK 


function [zros]=muller(fn,nprev,np,maxit,ep! ,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,ep | ,ep2,fnreal,zeros) 

% This function is called by the driver program mulrun.m 

% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, 1.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 

% by S. D. Conte and Carl de Boor, "Elementary Numerical 

% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

Jo 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
% no better estimates are available. 

% ep|: relative error tolerance on the root. 

% ep2: error tolerance on the function. 

Y% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% f{nreal:.a logical number; if TRUE, the program forces all approximations to 


% all the roots to be real. This makes it possible to use this routine 

% even if f(x) is defined only for real x. 

Jo Here we can take fnreal=1, 1f we only want to fine the real roots; we 
% can take fnreal=0 if we want to find both the real and complex roots. 


fn=’ string !4a’; 
eps |=max(ep1,1.e-12); eps2=max(ep2,!.e-20); 
for i=nprev+l:np; kount=0; 

while | 


% Compute the first three estimates for zero as 
% zero, zero+0.5, zero-0.5 


zero=zros(1i); h=0.00001+)*0.00001; zp=zero+h; zm=zero-h; 
[fzr,dvdflp,kount,zros]=dflate(fn,zp,1,kount,zros); 
[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 
hprev=-2.0*h; dvdfl p=(fzrprv-dvdf1 p)/hprev; 
[fzr,fzrdfi,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 


Yo Do while kount<maxit or h is relatively big or fzr=f(zero) 
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ae 


% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 


divdfl=(fzrdfl-fzrprv)/h; divdf2=(divdfl-dvdf1p)/(h+hprev); 
hprev=h; dvdflp=divdf1l; c=divdfl+h*divdf2; 
sqr=c*c-4.0*fardfl*divdf2; 
if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<O, den=c-sqr; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfi/den; fzrprv=fzrdfl; zero=zeroth; 
if kount>maxit, break, else 
while 3 
(fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if zros(ij==zero* 1.001, break; end 
To Check for convergence 
if abs(h)<eps1 *abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
Jo Check for divergence 
if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end; % end kount <= maxit 
if zros(i)==zero*l .001, break; end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
ZiOS\ I=Zero. 
end; % for loop 








Jo MULRUN 


J 2 2 2K 2 OG IS 2K OI 2 of IK OE 2k OE Of OK IK OIE IK OI IS IK OE OS OK OK IR IK OI OK HE 28 OIE 2S OK 2k 2 OK 


x=0; 
y=0; 
global R 
% input specs for Muller 
fOr N= bie 1 
R=500*(n-1); 
(Oe 
% tolerance, name of function that computes T12 etc: 
zrs = 104; % Starting guess for the root, can be real or complex 
fn=’stringl4a’;  % name of program that computes T12 
maxit = 500; % max # of iterations 
ep! "= le-5; Ytolerance on the root (start with 
%o values, le-2or so.) 
toliem = le-5: % relative tolerance on T12 
fnrea] = 0; % = 0 complex search, = 1, real search 
Metey — J leave these as 0 
io) = % and 1 
ep2 = tolfcn*string 14a(zrs); % This sets the tolerance on T12 to be 
% \le-> of T12 estimated at your 
Jo guess 


% Call Muller 


root = muller(fn,nprev,np,maxit,ep1,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 


subplot(3 11) 

plot(x,imag(y)) 

grid 

title’ mode, dampers at mode 3 antinodes’ ) 
xlabel(’R’) 

ylabel(’imag(nat freq)’) 

subplot(3 12) 

plot(x,real(y)) 
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grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 
subplot(313) 

plot(x,zeta) 

grid 

xlabel(’R’) 

ylabel(’ modal damping ratio’) 


9] 





Jo DFLATE 


% This is the subroutine needed by muller. 

function [fzero,fzrdfl,kount,zros ]=dflate(F,zero,i,kount,zros); 

% (*fzero,fzrdfl,kount,zeros* )=dflate(F,zero,i,kount,zeros) 

% This function which is called by function muller is translated 

% from the Fortran routine given by S. D. Conte and Car! de Boor, 
% "Elementary Numerical! Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 
% Engineering Mechanics Department, University of Alabama. 


kount=kount+1; fzero=feval(F,zero); fzrdfl=fzero; 
if i<2, return; end 
for j=2:1 
den=zero-zros(Jj-1); 
imabsiacni==0.0 
ZGOS(1)=Zeroa 100 1: remecn 
elise 
fzrdfl=fzrdfi/den; 
end 
end 





HANG 


The following is a typical sequence for using the model HANG to determine 


system natural frequencies and mode shapes and optimize system damping: 


1) In MATLAB, change directory to HANG. 
2) Edit STRINGI for desired absorber properties and STRING2 for desired cable 
properties. 
3) In STRING4 and STRINGI, ensure damping constant is set to zero. Run STRING4 
to plot t}2 of [Toveran] versus frequency. The zeroes of this plot are the undamped system 
natural frequencies. 
4) If desired, plot undamped system mode shapes as follows: 
a) For mode of interest, enter natural frequency in MATLAB window. 
b) Run STRING6 to plot mode shape. 
c) Repeat for other modes of interest. 
5) Perform absorber optimization for individual modes as follows: 
a) Edit MULRUN to iterate over desired range of damping constants. 
b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 
6) Set damping constant to global in STRINGI. Run MULRUN to plot real and 
imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant 1s that which maximizes damping ratio. 
7) Plot damped system mode shape as follows: 
a) Set damping constant to the optimum value in STRINGI. 
b) Enter the complex natural frequency in the MATLAB window. 
c) Run STRING6 to plot mode shape. 


8) Repeat steps 5 through 7 for other modes of interest. 





Jo STRINGI 


% This file contains the transfer matrix for a hanging mass spring 
% dashpot vibration absorber. 


%o Initial data 


IM=5>:- 
K=68.45; 
R15: 


t=]; 

UGle2 =O: 

t(2, L)=-(w°2*M*(i*w*R+K))/(K-w°2*M+i*w*R); 
(22 )=1- 
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%o STRING2 
% This file calculates the transfer matrix for a string of length L. 
JoInitial data 


ft —s000: 

L=28: 

p=1.0; 

c=sqrt(T/p); 

K=urc: 
t(1eip=cos(k*L); 

(alee =sin(k *Ly/{k* 1); 
ee) —-K* 1 =sin(k*L); 
122 =cos(k 1): 


oS 








Jo STRING4 


% This file plots the overall T(1,2) for a system with 3 components. 
W=0; 
Troe 2=0: 


==) 
1Omn—1- 200: 
Welln: 
W(n)=w; 
string2 
let 
ot 
string 1 
est: 
Tieet=—tel *T2*T3; 
heel 2(n)=Ttot( 2 
end 
plot(W,abs(Ttot12)) 
grid 
Xlabel( w(rad/sec)’) 
ylabel( T(1,2)’) 
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Yo STRING6 


%o THIS FILE PLOTS THE MODE SHAPES OF A STRING 
Jo WITH A MASS SPRING DASHPOT DAMPER AT ITS % 
BEFORE RUNNING, 

Jo ENTER FREQUENCY IN THE MATLAB WINDOW AND 
Jo ENTER DAMPING IN STRINGI. 

y=0; 

Oe 

i =5000: 

p=1.0; 

c=sqrt(T/p); 

k="v/C; 

x=022 1-50: 


fon n=12251 

UG = COStne” x4 i1))e 

eh? =sinke xumtk* 1); 
C1? een sintk= x(n); 
Ui. 2 J=COS(K™ x( 11); 
y(n)=tl(1,2); 

fit 292): 

end 


String] 

(Sl) )=cosi25-k): 
sale) Sim 2onk)/ (KW), 
iS(22)—COsi25 ah 
P52) =-k* cin 25°). 


for n=252:501 

ci) =COse(x(m)=25))): 
ily) —Ssimuieotx (ino) /(K* | ): 
teat =k esi k | x(m)-25)): 
(2.2 aC OS( men) -25)): 
M=tso tte 

y(n)=tf(1,2); 

fj=ti(2.2): 

end 


plot(eGh2s! reall. 251)),x( 2527501) real(y(252:501))) 
grid 


ai 


CENTER 





xlabel(’x location’) 
ylabel(’ real(displacement mode shape)’) 
title’ 3rd MODE MODE SHAPE’) 


to 
i 


9S 





‘ = 


“0 | | | | 


% STRING14 
Yo This file computes T(1,2) for a system with 3 components. 


function T12=string14(w) 
string? 

(lee 

T3=1; 

string | 

an 

Att le te 

JZ See 


oo 





an! 


Jo MULBER 


function [zros]=muller(fn,nprev,np,maxit,ep1,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,ep1,ep2,fnreal,zeros) 

% This function is called by the driver program mulrun.m 

% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, 1.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 

% by S. D. Conte and Carl de Boor, "Elementary Numerical 

% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

Jo 

% Yo find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

%e zros: an array of initial estimates of all desired roots, set to zero if 
Jo no better estimates are available. 

% ep: relative error tolerance on the root. 

% ep2: error tolerance on the function. 

% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% f{nreal: a logical number; if TRUE, the program forces all approximations to 


Jo all the roots to be real. This makes it possible to use this routine 

J even if f(x) 1s defined only for real x. 

Jo Here we can take fnreal=1, if we only want to fine the real roots; we 
Jo can take fnreal=0 if we want to find both the real and complex roots. 


fn=’string 14’; 
eps l=max(ep1,l.e-12); eps2=max(ep?2,1.e-20); 
for i=nprev+1:np; kount=0; 

while 1 


% Compute the first three estimates for zero as 
% zero, zerot+0.5, zero-0.5 


zero=zros(1); h=0.00001+j*0.00001; zp=zero+h; zm=zero-h; 
[fzr,dvdf1p,kount,zros]=dflate(fn,zp,i,kount,zros); 
[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 
hprev=-2.0*h; dvdflp=(fzrprv-dvdf1p)/hprev; 
[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 


Y Do while kount<maxit or h is relatively big or fzr=f(zero) 


100 


= 





ee 


% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 


divdfl=(fzrdfl-fzrprv)/h; divdf2=(divdfl-dvdf1 p)/(h+hprev); 
hprev=h; dvdfi p=divdf1; c=divdfi+h*divdf2; 
sqr=c*c-4.0* fzrdfl*divdf2; 
if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<0O, den=c-sqr; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfi/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 
(fzr,fzrdfi,kount,zros ]=dflate(fn,zero,i1,kount,zros); 
1zros()=—zero~ 1.001, breakvend 
%o Check for convergence 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
% Check for divergence 
if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end: Yo end kount <= maxit 
if Z50s())——zero 1-001) break: end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<eps1*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while | 
ZIOS() =ZEro. 
end; % for loop 


10} 





Jo MULRUN 


Oy 5 RR RR ROR AOI CR IO CK 


x=0; 
y=, 
global R 
% input specs for Muller 
eee 
= Sia 
x(n)=R; 
% tolerance, name of function that computes T12 etc: 
zrs = 4.07; % Starting guess for the root, can be real or complex 
fn = ’string14’; % name of program that computes T12 
maxit = 500: % max # of iterations 
epl =  ile-5; Ytolerance on the root (Start with 
% values, le-2or so.) 
tolfcn = le-5; % relative tolerance on T12 
fnreal = 0; % = Q complex search, = 1, real search 
fire ve—U, % leave these as 0 
np = 1; % ayikal II 


ep2 = tolfcn*string14(zrs); % This sets the tolerance on T12 to be 
J le-5 of T12 estimated at your 
Yo guess 


% Call Muller 


root = muller(fn,nprev,np,maxit,ep1,ep2,fnreal,zrs); 
y(n)=root; 
end 
zeta=imag(y)./real(y); 


subplot(311) 
plot(x,imag(y)) 

grid 

title? 4th mode’) 
xlabel(’R’) 

ylabel( imag(nat freq)’) 
subplot(312) 
plot(x,real(y)) 
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lower 


Starting 








grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 
subplot(313) 

plot(x,zeta) 

grid 

xlabel(’R’) 

ylabel(’ modal damping ratio’) 
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Jo DEFLATE 


% This is the subroutine needed by muller. 

function [fzero,fzrdfi,.kount,zros]=dflate(F,zero,i,kount,zros); 

% (*fzero,fzrdfi,kount,zeros* )=dflate(F,zero,i,kount,zeros) 

% This function which is called by function muller is translated 

% from the Fortran routine given by S. D. Conte and Carl de Boor, 

% "Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 


kount=kount+1; fzero=feval(F,zero); fzrdfl=fzero; 
if 1<2, return; end 
for J=2:1 
den=zero-zros(j-1); 
if abs(den)==0.0 
zros(1)=zero* 1.001; return 
else 
fzrdfi=fzrdfi/den; 
end 
end 
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Appendix B: Transfer Matrix Library 


The following appendix is a direct reproduction of the transfer matrix library 
derived by Li Li [5]. It is included here so that this thesis can serve as a stand alone 
reference on transfer matrix analysis 
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Appendix A 


Transfer Matrix Library 


Definitions of the State Vectors and Transfer Matrices: (see Section 2.1) 


A.1 Strings 


1) A string with varying tension 
The WKB solution to a string with varying-tension (see Fig. A-1) is given in [20] 


as 


y Ciel Se Cen A. 
I= or pT at. (A.1) 


The transfer matrix can be found from this solution as 











t t 
r= 11 %i2 (A.2) 
to1 te2 
where 
tui= (42 )tcosW + Tict gin W 
(pi teat 2a) 
t 1 
12 We psinW 
ae 17. : To ee : 
toy = [SE( gE )t — 3(41)s]cosW — [w(p?T:-1T:)* + ge) sinW 
: (e2T,_17,)4 
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i-1 


2) 


Figure A-1: A string with varying tension 


ye ae ee. zy 
22 = (7) cosW EEC 


where W = w fi Wee If the tension varies linearly, then W = Pe /T; = 


Vian Heat =n", 


2) A string with constant tension 
For a string with constant tension (see Fig. A-2), with k = vj eee eee the 


transfer matrix is 


ees coskl, pz sinkl, (A.3) 
—kTsinkl,  coskl. 


where [, is the distance between the two state vectors at 2 — 1 and 1. 
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Figure A-2: A string with constant tension 





Figure A-3: A general, rigid lump 


A.2 Rigid Lumps 


1) A general, rigid lump 
Section 2.5 gives the detailed derivation of the transfer matrix for a general, rigid 


lump (see Fig. A-3), the result is 


1 Ji Gis =- l) + T._4il-— Mw? |r? + (L. = DP] lig 


~Muw?(T(L.—1 4 T.yl-— Mw?r*) (LL. -1) + T_yl — Mw? (r? + 2] 
(A.4) 


F; = 


>| 


where A = T,(L, — 1) + T;_,l — Mw?[r? — (2, — Dj. 
2) Application 1: a uniform bar 
For a rigid, uniform bar (see Fig. A-4) with length J, and linear density p;, then 
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p 
il l 
ta <— 
-~—» 


Figure A-4: A rigid, uniform bar 





Figure A-5: A uniform, solid sphere 


a 7 M = pil,, and r* 


[? oe 
EE The transfer matrix is 


4(3 — F a 
i oll (3 — Jy) T (A.5) 
26+) | —mw(12-) 4(3- 4) 


where J, = 4" and T is the average tension if there is a difference in tension at the 
two ends of the lump. 


3) Application 2: a solid sphere 


For a uniform, solid sphere (see Fig. A-5) with diameter D and volume density 


3 ° e 
Do uleneL. = Hern — as HV eens and r? = De The transfer matrix is 


| Bias } 7 (A.6) 
20+ 31, | _Mu?(20-2/,) 20-7, 
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Figure A-6: A cone 


7) - e e @ ® « e 
where J, = wea and JT is the average tension if there is a difference in tension at 


two ends of the lump. 


4) Application 3: a cone 


For a cone (see Fig. A-6) with height h, bottom radius r,, and volume density p,, 
then L. = i). l= Sh, \Va— sg uLeD Reed se 3(r24+2h7) 


= ee The transfer matrix is 


7 1 L—Ie(r§ + $A?) rata . 
eS ee nT (A.7) 
L— ders — 2°) | _ ages? — Fe(r? + 1h?) 1 — I(r? + 4h?) 
w? 
Gene Js = gaa: 


5) Application 4: a point mass 


The transfer matrix for a point mass with mass M can be obtained directly from 


the transfer matrix for the general, rigid lump with D,. =l=r=0: 
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c=l],+15 d=I3+l4 e=Iel¢ 


Figure A-7: A general, rigid lump with damping devices (type I) 
A.3 Damping Devices 


For a rigid lump with a damping device, the transfer matrix is defined the same 


way as that for a rigid lump without a damping device, 


y E C11 by2 y (A.9) 


ey C21 t22 ey, 


t t—1 


The expressions for the elements of the transfer matrix are as follows: 


ae Fa 3 a a5 
Pee ee ee ee eM 4 IS 2s 
l — ¥3°¥6 1 — 73776 1 — 736 ee ee 
(A.10) 
For the system shown in Fig. A-7, y’s are given as follows: 
_ ga(91912 — 93911) — 95(92911 + 9199) (A.11) 
912(9296 — 9498) + 99(949s + 9396) 
= 95(99 — 92910) — ga(93910 + 912) (A.12) 


oo 
912(9296 — 9498) + 99(949s + 9396) 


11] 





___ ga( 9397 — 9s) + 912(g297 + gs) 
— (95) F 9129297 F 98) (A.13) 
912(9296 — 9498) + 9o(949s + 9396) 


= 91(95918 — 98913) + 917(929s + 9398) (A.14) 
gis( 9397 — 9s) + 913(9297 +98) + gi4( 929s + 9398) 


98913 — 915(929s + 9398) — 95918 
“5 = ) (A.15) 
gis(9397 — 9s) + 913(9297 + 98) + g14(929s + 9398) 


et 913( 9296 a 949s) = 916( 929s a5 9398) = 918( 949s a 9396) (A 16) 
gis(9397 — 9s) + 913(9297 + 98) + 914(929s + 9398) 


where g’s are given by 


y= w?M,(1 = 1), g2 = w?( Myo ~ M,‘+), 93 = w?( M4 + M,*+), 
gg = w?M3(1— 8), gp = Btw? Mae +B gg = 4% + w? Meals (42 — 1), 
Gi. -8, 9s = 2 go = Btw? ik +o. 


Jio = ¢, 911 = a a w? My12(1 = .), G20 a 


a= = = w? My ‘4s “Le Pa 914 = la, Qis = Is, 

ore — w? M3l4(1 = ‘4 ) se oe nNi= on i w? M,13(1 a iL), S18. =a a w? My oS 5° o. 
a, = Til, + Ail, + Kya? —w7], + wR, 3} 

a. = H,ly+ Ayl3 + Kya? + Ka? + ju( Rb? + R2b3) — wl, 

a3 = Tile + H,ls + K2a3 — w?I3 + jwR2b3 


B, = Kya? + jwR,b?, B, = Kya} + jwR2d3 


where H; and H, are the forces along the structural length at the left and right 
joint, respectively, which can be found by statics. 


For the system shown in Fig. A-8, +y’s are given as follows: 


_ —(M,l+ M2f)a, — M,ca3 + w?M, M2 fe? —M,fh? — M,ld* + (a,l + a3h)/w’ 


Mia, + (Myd-+ Mzh)ag +o? M, Maced’ =" Mea, + (Mid + Myh)as + 0? My Mae 
(A.17) 








Vg.n = GaGis + 93913 — 915(9293 + 9393) — 95918 — 959\8) 

Yen = 913(9296 — 9498) + 913(9298 — 9498) — 9ie(G29s + 9398) — 916(9295 + 9398) 
+93a(949s + 9396) + 918(9495 + 9396) 

Vad = 918(9397 — 9s) — 91895 + 913(9297 + 9s) + 91393 + 914(9295 + 9395) 

15,d = Yéd = pane 

gg = jwb?(2 + 4), 96 = jwb?=, 9 = jwd? 5, 

gg = jwh?(2 + 5), 91 = Iw 2, ie = Jw? 5, 

913 = jwb? (4 ag *); Jig = jwb?*, 97 = jude, 

Jig = jub?(F ate A) 


For the system shown in Fig. A-8, y’s are given by 


 _ jw? {{M,(l+e)+Me f}ry at[Mi d+Me2(eth)]yi 1} ~! oe ju? {— Py. at(MidtM2(e+h)] 72,2} 
a 


N Vd 13Wd 


jw? { SPs 4+{Mid+M3(e+h)]13,n} = jw? {—[My (I4e)+Mz f)y4.a- Sy 74,0} 


SS eee = 
%3 13 4 4 Vaid 


jus ub? (2 (ys.a—e.n)) fwd? {- [My d+ Mo (e+h)]¥6,4-*F 76,0} 


» Ye = Ted 





aren 


A.4 Supports 


For the support shown in Fig. A-9, the transfer matrix is 


_ L.-2z 0 
F.= 5 z+0,L A.26 
Tow? —T;(Lewz)-Tinsz ez | -_— ( ) 
(Le—-z)z L.-z 


where J, is the mass moment of inertia of the rigid lump about the pinned point. 


For the support shown in Fig. A-10, the transfer matrix is 





_ L.-2z 0 
i . eee Oe (A.27) 
Iow? -T(L.-2) z 
(Le—-z)z ae 


where J, is the mass moment of inertia of the rigid lump about the pinned point. 





ie aa Z — hinge —> 7. 


Figure A-9: An intermediate support 


Figure A-10: An end support 
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Appendix C: 


Analysis of a Single Degree of Freedom System Using the 
Transfer Matrix and Complex Natural Frequency Approach 


The Svstem: 


— — oo vo 









! K (N/m) 
! F; 
Va ! M (kg 
R (N/(m/s)) z 
omer Es ! 
— ae ! 
The Transfer Matrix Relates the State Vectors at Points | and 2: 
lie Ti Tiz | x: 
F:| | Tx To: |! Fi (1) 
Derive the Transfer Matrix of the Svstem Between Points | and 2: 
-F, = K(x2 = X1) + joR(x2 - X1) a (K + JOR)(x2 = X1) 
solving for x2: 
xooa -F,/(K + J@R) ao) (2) 


F, + F>= M(d°x2/dt*) 

F, + F) = -0’Mx> 

writing F> in terms of F; and x;: 

F, = -F, - @°M(-Fi/(K + joR) + x1) 


F2 = -@° Mx; + (o M/(K + joR) -1)F, (3) 


Equations 2 and 3 result in the transfer matrix: 


Ta Tz Se -1+@°M/(K+joR) 


7 i, . -1/(K+jaR) 


Applying Boundary Conditions: 


X; =0 due to system being attached to wall 
F,=0 because we want to study the free response of the system 


so, equation | becomes: 
X2 Tu Tiz | 0 
obs tale 
ylelds the characteristic equation: 
Ta =0 
-1+@°M/(K+joR) =0 
-0o°M + joR+K=0 (4) 


The Complex Frequency Root: 


applying quadratic formula with 4KM < R’: 


V4KM-R°> oR 


= ae laa 
2M 2M 


= Oreal + JOunasiee (5) 
From Section 3.3.1 of thesis, a damped harmonic response can be written: 
x(t) — 7 COS(Mreat t) pee 


V4KM -R°- R 
peers || CXD on (6) 


x(th=A 
x(t) co OM 
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Comparison With Traditional Method: 


In the traditional method, the characteristic equation is the same as 
equation 4: 
-°M + joR+K=0 


but the complex root is not solved for. Instead the following are defined: 


|K 
Oo= VF the undamped natural frequency (7) 





R 
¢ = 5 az the damping ratio (3) 
Md 


Od = Wo eee , the damped natural frequency (9) 

The displacement solution is then written as: 

x(t) = A cos( @at ) exp(-¢ @ot) (10) 

When equations 7, 8 and 9 are used to write equation 10 in terms of 


K, Rand M. The result ts: 


V4KM - R- R 
ener ware CN) rye i) 


X(th=A co OM 


which is the same as equation 6, the result obtained using the complex 


frequency approach. 
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